clear;

delta_exac = 0;
delta_predict = 0;
cnt = 0;
filename = "case9";


mpc = loadcase(filename);
mpc = ext2int(mpc);
mpc.bus(:,3) = 0;
mpc.gen(:,2) = 0;

gen_buses = find(mpc.bus(:, 2) == 2 | mpc.bus(:, 2) == 3);
load_buses = find(mpc.bus(:, 2) == 1);
%randomize the raw data
gen_selected = gen_buses(randperm(length(gen_buses), round(0.3 * length(gen_buses))));
gen_selected_idx = find(ismember( mpc.gen(:,1), gen_selected));
load_selected = load_buses(randperm(length(load_buses), round(0.3 * length(load_buses))));

a = normrnd(0, 0.5);
mpc.bus(load_selected, 3) = (1+a) * mpc.bus(load_selected, 3);
mpc.bus(load_selected, 4) = (1+a) * mpc.bus(load_selected, 4);
b = normrnd(0, 0.3);
mpc.gen(gen_selected_idx,6) = (1+b) * mpc.gen(gen_selected_idx,6);

%predict whehter the system is stable
Ybus = full(makeYbus(mpc));
Ybus = imag(Ybus);

BGG = Ybus(gen_buses, gen_buses);
BLL = Ybus(load_buses, load_buses);
BLG = Ybus(load_buses, gen_buses);
[~, genidx] = ismember(gen_buses, mpc.gen(:,1));
VG =  mpc.gen(genidx,6);
V_L_star = -inv(BLL) * BLG * VG;
Q_crit = 0.25 * diag(V_L_star) * BLL * diag(V_L_star);
QL = mpc.bus(load_buses, 4)/mpc.baseMVA;
temp = Q_crit\QL;
Delta = norm(temp, "inf");

if Delta <= 1
    suc = 1;
else
    suc = 0;
end

%exact solution
mpopt = mpoption('out.all', 0, 'verbose', 0);
[result, success] = runpf(mpc, mpopt);

if success == 0
    disp("fk");
else
    bus_voltages = result.bus(load_buses, 8);
end

%whether the prediction is right or not
if success == suc
    cnt = cnt + 1;
    delta_predict = delta_predict + (1 - sqrt(1 - Delta))/2;
    delta_exac = delta_exac + norm(abs(bus_voltages - V_L_star)./V_L_star,"inf");
else
    disp("It fails.")
end

accuracy = cnt / 10;
delta_exac = delta_exac / cnt;
delta_predict = delta_predict / cnt;
condition_accuracy = abs(delta_exac - delta_predict)/delta_exac;


accuracy
delta_predict
delta_exac
condition_accuracy

